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A new technique, root projection (RP), 
is given for quantitative deconvolution of 
causal time series in the presence of 
moderate amounts of noise. Deconvolu- 
tion is treated as a well-conditioned but 
underdetermined problem and a priori 
information is employed to obtain com- 
parable noise reduction to that achieved 
by singular value decomposition (SVD) 
techniques while providing more accu- 
rate frequency information about the in- 
verse. Two detailed examples arc given. 
The first gives noise analysis for alter- 
nate methods for deconvolution with a 
Gaussian kernel. The second example 



presents a model acoustic emission 
transducer calibration problem with typi- 
cal noisy and incomplete output data. 
This example is treated by the use of a 
robust cross-cutting algorithm combining 
both the RP and SVD methods. 

Key words: acoustic emission; causal 
time series; cross-cut algorithm; decon- 
volution; FFT deconvolution; Gaussian 
deconvolution; root projection; signal 
processing; singular value decomposi- 
tion; time series; transducer calibration. 

Accepted: Februaiy 25, 1991 



1. Introduction 

Because of the almost ubiquitous presence of 
Green's functions in linear physical theories, de- 
convolution (or inverse filtering) is the most com- 
monly encountered special inverse problem 
associated with many characterization problems in 
NDE, both electromagnetic and elastic. A difficulty 
with inverse problems is their susceptibility to ill- 
conditioning; and, as is well-known, deconvolution 
problems can range in difficulty from those simply 
solvable by forward substitution to those for which 
no known method holds any hope of solution. 

Some of the most powerful techniques that can 
be applied to deconvolution are those, such as reg- 
ularization and singular valued decomposition 
(SVD), arising from general linear inversion meth- 
ods. However, there are other methods available 
which take advantage of the special properties of 
the convolution process. We shall review briefly 



some of these methods, compare them with an 
SVD method especially adapted for deconvolution 
and present a new approach to deconvolution rest- 
ing on the concept of root projection which, when 
combined with SVD, offers the possibility of ap- 
proaching currently intractable deconvolution 
problems. 

Two examples— both ill-conditioned— will be 
presented, and signal-to-noise concepts will be in- 
troduced to sharpen the SVD and root projection 
deconvolution (RPD) methods. These will then be 
compared with more standard deconvolution meth- 
ods. The first example involves convolution with a 
Gaussian kernel in which the entire output is given, 
but corrupted by varying amounts of noise. In the 
second example the convolution kernel is derived 
from the actual response of a "standard" acoustic 
emission transducer in a calibration experiment. As 
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is typical in many such experiments, the entire out- 
put is not available. So one is faced with deconvolv- 
ing a one-sided Icernel from an incomplete 
one-sided time series output. To treat this prob- 
lem, the SVD and RPD methods are combined 
into a more robust cross-cutting algorithm which 
provides both superior time and frequency infor- 
mation about the inverse. 



2. Deconvolution Methods 

Full convolution of two finite-length, real-valued, 
one-sided time series may be expressed as 



a:N*b:M = (a *b ):(M +N), 



(1) 



where a is a causal time series with elements 
ao,...,aN, 6 is a causal time series with elements 
bo,...,bM, and where a *6 is a causal time series with 
elements' 



(a *b )„ = XkUk b„.k,n= 0,-,M +N. 



(2) 



The fact that a:b or b:M are truncated at finite 
length makes no difference to the first min(M,N) 
terms in the output. 

In the deconvolution problem one is given the 
series a:N and some part of the right hand series in 
eq (1), usually M terms, which is possibly corrupted 
by noise. For most applications one assumes that 
M=N, if need be by zero padding the shorter trun- 
cated series, and one writes the deconvolution 
problem in the form 



a:N*b:N = v:N 



(3) 



where a:N and v:N are assumed known and p:N is 
to be reconstructed as closely as possible. The se- 
ries a:N is called the kernel series. In Y transform 
terminology eq (1) may be paraphrased as 



a(y)P(y) = v(y)+y''^'u(y), 



(4) 



where u:(N — l) is usually an unknown time series. 
In the notation used in eq (1) with N <M, if a:N 
is thought of as an A^th order difference operator, 
then eq (1) becomes a high order difference equa- 
tion with the first K terms associated with the 
boundary conditions at the causal boundary and 
the last N—K — 1 terms associated with the 



' We use here the modified Dirac bracket and y transform nota- 
tions employed in references [1] and [2]. 



boundary conditions at the transient boundary. 
Treating eq (3), then, has the interpretation of 
solving a one point boundary value problem, while 
eq (4) may be thought of as solving the two point 
boundary value problem where the influence of the 
transient boundary conditions on the solution may 
be studied and these conditions possibly restricted 
using a priori information. 

In Y transform terms the inversion of eq (4) 
merely amounts to finding the coefficients of the 

Taylor's series of ivfy) +u(y)y'^*^ j/a (y) about the 

origin, and the first N + 1 of these coefficients are 
the same as those of v(y)/a(y) so long as flo^^O. 
From this point of view the expansion of v(y)/a(y) 
about the origin is only valid within the radius of 
convergence of the Taylor's expansion. Since v(y) 
is a truncated series and has noise, its zeroes will 
not include those of a(y), so that this expansion is 
generally only valid to the innermost zero ofafy). 
Let the modulus of the smallest zero of a (y) be 
/■fl. If we use a circle of radius r<rQ, then we can 
write the Cauchy formula for the Taylor's expan- 
sion coefficients as 






(5) 



Approximating this integral by numerical evalua- 
tion on r times the Mth roots of unity gives: 






2irW\ 2irfa' 



(6) 



This formula is exact when O^n^N <M if 
v:N = (a*j3):N. Otherwise, all coefficients /3„', 
n' =m (mod M) of the Laurent expansion of v(y)/ 
a(y) will be aliased onto /3„. 

If one defines the "radiused" time series 'ar:N, 
where (ar)„ =a„r", n =0,—,N, then eq (6) is merely 
a DFT inversion formula on the quotient of the 
DFT's of the radiused time series vr:N by 'ar:N 
where the n th term of the inverted series is multi- 
plied by the expansion factor of r "". We refer to 
this technique implemented with the FFT al- 
gorithm as radiused FFT inversion. The factor r~" 
in eq (6) expands noise exponentially rendering 
this algorithm generally useless in its raw form. The 
algorithm can, however, be employed in conjunc- 
tion with root projection. That application will be 
discussed in detail shortly. 

If one attempts to employ eq (6) for r<roto re- 
duce the noise buildup, one obtains the expansion 
of the wrong Laurent series form of v/a, a form 
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yielding a non-causal time series whose terms de- 
pend on the residues of u(y) at the roots of a(y) 
inside r. In moving r from the inside to the outside 
of any one of these roots one changes the contribu- 
tion due to that root from a causal time series with 
a factor of the residue times r~" into an anticausal 
time series with terms of the type residue times 
r"~K In a DFT context these terms also wrap 
around through negative times to produce time 
aliasing, especially for large positive times. The ra- 
pidity of convergence of the partial fraction series 
associated with a given root ya depends on that ra- 
tio yjr or rfya which is less than one.^ 

In matrix notation the deconvolution problem is 
almost exclusively formulated in the parallel form 
of eq (3): 
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while eq (4) takes the form: 
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(8) 



Equation (7), while completely determined, is gen- 
erally ill-conditioned. Equation (8), on the other 
hand, while underdetermined, is generally well- 
conditioned and offers new algorithms for dealing 
with the deconvolution problem. 



^ In the ideal continuous case there is no truncation or aliasing 
problem associated with v(y), only the noise problem. In that 
case the Fourier transform for p(y) should be entire in the up- 
per half plane, and the exponential noise buildup and non- 
causal representation arguments against the continuous Stokes' 
deconvolution formula for dividing the appropriate Fourier 
transforms are still valid [3]. 



The simplest method for solving eq (7) is by for- 
ward substitution, which is equivalent to formal di- 
vision of the y transforms v(y)/a(y) starting at the 
constant term. This is generally the most unstable of 
inversion methods. The matrix equivalent of zero- 
padded and radiused FFT inversion is obtained by 
considering inverse-radiused circulant matrices. 

From a real-valued time series c;/Vwe build the 
(/V + 1) X (/V + 1) inverse-radiused circulant matrix 
Cr defined by: 



(Cr)u =r C(mo6{k-l,N+l)),k,l = 0,"*,/V. 



(9) 



It is then a simple matter to show that the vectors 



r 'e 



—l-nmi 
—4Tmu 



N-¥l 






(10) 



are eigenvectors for C. That is: 

-iTTmi r / Itrmi \T "1 

[a]r-^e^^:/V=[(e;^:/v) c:N\ 

(-2imi \ 



r 'e 
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(11) 



where c:N is the /V-f 1 point DFT of c:N. In addi- 
tion, the radiused DFT provides a decomposition of 
the discrete S function: 

re^;/Vjr-^e7^:yV=(/V+l)S„„. (12) 

Then, the inversion of the equation: 
CrP:N = v:N (13) 

is easily carried out by eigenvector decomposition: 



/3;A/ = 






(14) 



To relate c:N to a:N, we set Ck =aij-'', k = 0,...,N, 
Ck = 0,k =N + 1,...,/V, which is the radiused form of 
the zero-padded kernel series. Then the inverse-ra- 
diused matrix Q has the form shown in eq (15) of a 
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(15) 



zero-padded convolution matrix with an extra ex- 
ponentially damped upper right hand Toeplitz cor- 
ner. Thus radiused FFT inversion solves the 
"damped" circulant rather than the convolution 
matrix problem. As seen from eqs (12) and (14), 
the answer obtained from this method will be cor- 
rect for the convolution problem if r^**'^ times the 
last N terms of the answer are effectively zero. 
That happy circumstance is usually frustrated by 
the exponential buildup in the eigenvectors of eq 
(10) together with the presence of noise compo- 
nents in v:N, as discussed previously in the Y trans- 
form context. However, referring to that 
discussion, one can conclude that if r < 1 and if 
there are no zeroes of a(y) with modulus less than 
r, then with sufficient zero padding FFT inversion 
will approach convolution inversion. 

Least squares techniques are the standard way to 
deal with noise problems. Frequently they are 
formulated in terms of minimizing the norm 
||i;;iV-A P:N\\. However, since v:N contains signal- 
plus-noise, applying this condition too stringently 
forces induced domain noise into /3;iV. What is re- 
quired is a technique to separate the range signal 
from the range noise either by applying a priori in- 
formation about the signal in the range or a priori 
information about the signal in the domain. In a 
simple sense the example given in reference [2] giv- 
ing rise to figure 4 accomplished this in the range 
using root projection when the complete informa- 
tion was available in the domain. We shall extend 
that concept below. 

Another way to deal with this difficulty is 
through regularization, in which the norm is modi- 
fied by adding another norm— usually a scalar 
times a domain norm— onto the above range norm 



to adjust the range residuals to be of the order of 
magnitude expected from the noise. The question 
is: Which norm should one add? The techniques 
for choosing an appropriate norm and carrying out 
this process efficiently lie somewhat outside the 
framework of this work, so that we shall not 
provide detailed comparisons with regularization 
methods here. 

For eq (7) the least-squares method is usually 
implemented by multiplying the equation by 
A^=^(A^' to give the equation 



A^A P:N =A^ v:N 



(16) 



where A'^A is a non-negative Hermitian matrix. 
The multiplication accomplishes two results: /) It 
forces the output vector into the range of A, 
thereby eliminating some of the components of the 
noise, and ii) It modifies the left hand matrix to a 
Hermitian form allowing eigenvector decomposi- 
tion methods and algorithms to be employed. How- 
ever, this approach suffers from two drawbacks: /) 
If^ has a simple structure— e.g., y4 is a convolution 
operator— efficient algorithms which take advan- 
tage of this structure may be lost, and ii) The ei- 
genvalues of A'^^A are effectively the squares of 
those of A (in a meaning to be described directly), 
thus greatly increasing the ill-conditioned nature of 
the problem. 

A closely related approach to the above is the 
singular value decomposition technique (SVD). In 
the SVD method one need not multiply the right 
hand side of eq (7) by A *. Rather, one represents 
A in the form 



A=U^DV 



(17) 
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where U and V are unitary matrices and £) is a 
diagonal matrix whose diagonal elements are called 
the singular values of A . The SVD method, while 
usually slower to compute than the method arising 
from eq (15) has the advantage of finding diagonal 
elements which are the square roots of those stem- 
ming from A'^A. This permits greater dynamic 
range for treating very ill-conditioned problems 
such as frequently occur when solving eq (7). For 
the deconvolution problem the SVD algorithm has 
been implemented using an efficient form of an- 
other algorithm due to Lanczos (in addition to the 
SVD decomposition) to permit a small dimensional 
approximation to the full SVD approach [4]. How- 
ever, this Lanczos/SVD approach, while more effi- 
cient may not permit accurate inversion especially 
when N is much greater than the dimension of the 
Lanczos subspace. 

We turn now to least squares algorithms associ- 
ated with eq (8). Calling ^2 the {2N + lxN + 1) 
matrix in that equation, one can say that A 2 de- 
scribes an embedding of the generally iV + 1 dimen- 
sional range of^ in 2N + \ dimensional space. In 
terms of the Euclidean inner products for both 
6^""' and ^'^^^ given by {u,v)=u'^v*, the adjoint 
of ^2, ^2^ isA2''=A2''=A2^, since^2 is real [5]. 
Given a vector c:2Nv/g can use A2^ to construct a 
least squares projection to map c:2N into the range 
ofA2. Once in the range, we shall see that inver- 
sion is well-conditioned and simple. 

The necessary and sufficient condition that c:2N 
lie in the range of ^ 2 is that c:2N lie in the orthog- 
onal complement of the null space of ^2^ [5]. That 



is: 



UA2H:2N = 0, then (p:2NYc:2N = Q. 



(18) 



To find a basis for the null space of A 2^, we con- 
sider the geometric sequence vector "J^ where y is 
any complex number. These vectors are eigenvec- 
tors of ^2^ satisfying 



A 2+ y::2N = (y::2NfA 2=a(y 'f::2N. 



(19) 



Thus, if a (y ) has distinct roots, we can identify the 
N dimensional null space of ^2^ as spanned by 



y^, where a (y,,) = 0, a = l,*",iV. 



The condition 17 then becomes the relation thatj'a 
are among the roots of c(y), or in algebraic terms 
a(y) divides c(y). From this point of view the exis- 
tence of CBR implies a set of "continuity" equa- 
tions that must be satisfied by time series in the 
range of a(y) under polynomial multiplication. 



The root projection discussed in reference [2] 
provides us with the required mapping into the 
range of ^ 2 [2]. Also, since there must exist a poly- 
nomial b(y) of degree N such that a(y)b(y)=c(y), 
we can find this polynomial simply by radiused 
FFT division with at least 2N + 1 points using any 
numerically reasonable radius near one. 

Thus we can formulate the root projection de- 
convolution (RPD) algorithm as an alternative 
least squares approach to deconvolution: 

i) Embedding kN in 2N + 1 dimensional space 

This is usually done by replacing an a priori 
range estimate with more accurate a posteriori 
range data where available. It can be done from 
the domain by convolving an a priori domain esti- 
mate with the kernel to give the a priori range esti- 
mate. Since the root projection of an a priori range 
estimate in this case is equal to itself, an equivalent 
method here is to subtract the a priori range esti- 
mate from the a posteriori range data, where avail- 
able, and fill out the rest of the range with zeroes 
to produce a reduced problem with a zero a priori 
estimate. The a priori domain estimate can then be 
added back on to the inversion estimate from the 
reduced problem to give an upgraded inversion es- 
timate. This latter approach can also be directly 
applied to introduce a priori information into the 
Lanczos/SVD method. 

ii) Applying root projection 

Here time weighted inversion is employed using 
weightings on the a priori and a posteriori parts of 
the range data reflecting the relative uncertainties 
of the a priori versus a posteriori parts of the data. A 
very high weighting factor applied to the a priori 
part of the data will cause the root projection to 
return a result with very small residuals in the a 
posteriori part of the range data. This means that 
the inverse fits both the range signal and the super- 
imposed noise more closely as one increases the 
ratio of weights of a priori to a posteriori points in 
the time weight file. 



(20) iii) Inverting using FFT division 



The entire FFT division operation is carried out 
on an A/ > 2N^ + 1 dimensional space. The basis for 
this space is the set of DFT vectors built from the 
(M - l)st roots of unify. However, only N of the M 
FFT components of any domain vector are inde- 
pendent, since there are M — (N + 1) conditions re- 



349 



Volume 96, Number 3, May-June 1991 

Journal of Research of the National Institute of Standards and Technology 



quiring all components after iV + 1 to be zero. In 
the range there are still only M - (N + 1) FFT com- 
ponents independent. This time there are N condi- 
tions requiring orthogonality to the geometric root 
vectors of the kernel and M-(N + 1) conditions 
requiring all image space components (2N + 1) to 
be zero. Thus, the basic range projection consists 
of projecting through the set of geometric root vec- 
tors supplemented by any simple basis of the 
M — (N + 1) dimensional space of vectors whose 
first 2N + 1 components are zero using the time (or 
frequency) weighted norm. The basic domain pro- 
jection is through the M-(N + 1) dimensional 
space whose first N + 1 components are zero. An 
independent orthogonal base for the range space is 
never built in root projection. 

Few practical deconvolution problems can end at 
this point. Because of the presence of range noise, 
any inversion process will induce a domain noise 
which frequently obscures some or all of the fea- 
tures of the inverse. The critical problem when car- 
rying out deconvolution consists of separating 
those aspects of the inverse which are "unequivo- 
cally" determined by the inversion process from 
those which are not. From the SVD or RPD point 
of view this means finding the signal-to-noise ratios 
in a series of orthogonal channels, accepting those 
aspects of the signal in the channels where the ra- 
tio is high and trying to supply a priori information 
in those where it is not. The SVD method provides 
an excellent means for optimal signal-to-noise sep- 
aration. RPD, which works with frequency infor- 
mation, often provides an excellent representation 
for testing and including a priori information. One 
can add or remove various hypothetical features 
from an a priori estimate and test the consequence 
of these changes. Since with adequate weighting 
root projection always forces the range residuals 
back to within the noise limits, the result obtained 
will be a possible solution. However, essential solu- 
tion features will be added and inessential ones 
taken away will not return. 

Both RPD and SVD offer the possibility for de- 
termining the relative probability of a particular 
deconvolution estimate. However, we shall not dis- 
cuss statistical algorithms in any systematic way 
here. Rather, we shall present two numerical ex- 
amples principally to compare the accuracy and 
flexibility of the SVD, RPD, and radiused FFT di- 
vision algorithms as inversion techniques and to 
present a cross-cut method for deconvolution 
combining Lanczos/SVD inversion (and filtering) 
with RPD and frequency space filtering. 



In the first case discussed the noise we add is 
produced by rounding off the range data to a cer- 
tain degree of accuracy, thus simulating an A/D 
conversion process. This noise while crudely uni- 
form in frequency amplitude is not independent of 
the output, since one expects a negative correlation 
with the output derivative. In the second case the 
noise added was produced by a random number 
generator. Error will be measured by using a rela- 
tive standard deviation. The ordinary standard de- 
viation of one file from another will be divided by 
the generalized geometric mean of the non-zero 
terms of the reference file to give an order of mag- 
nitude correction to the error without being overly 
sensitive to "spikiness" in the reference file. Errors 
will only be compared for the same reference file. 



3. Examples 
3.1 Example I 

We return to the example arising from the con- 
volution of the 101 point Gaussian of figure 1 in 
reference [2] with the first 700 points of the experi- 
mental waveform of figure 2.4a in reference [1]. 
This example has the form of eq (1) with full range 
information; only noise has been added. This ex- 
ample illustrates one way of looking at M +N di- 
mensional deconvolution in terms of root filtering. 
The 8 bit noise added on in that example was an 
"arbitrary" element of M+N + 1 dimensional 
space. 

The projection condition placing the signal in the 
range of the kernel reduced the "power" of the 
noise by one eighth in this case. Another 700 de- 
grees of freedom in the noise are connected to the 
roots of the Y transform of the answer, )3(y). Al- 
though usually impractical, each piece of a priori 
information describing a root of )3(y) removes one 
degree of freedom from the noise. When the roots 
of both the answer and the kernel are known, only 
the magnitude of the answer is undetermined. This 
is reflected in the fact that the only degree of free- 
dom left for the noise is associated with the one-di- 
mensional space generated by the convolution of 
the kernel and the answer, any element of which is 
orthogonal to the geometric series root vectors 
from both the kernel and the answer. 

Two levels of range noise will be imposed for this 
example. We first consider the case where only 24 
bit noise is added to the output, thus reducing the 
output to "single precision" accuracy. 

Even with 24 bit accuracy the use of the forward 
substitution (or "real time deconvolution") method 
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for this problem produces errors in excess of 400% 
by term 100. The errors increase exponentially af- 
ter that point. 

Because most of the roots of the Gaussian lie on 
the unit circle, we carry out our analysis on a circle 
inside the unit circle. We choose a radius of 0.997 as 
reasonable for dealing with a series of length 700 
and a zero-padded FFT of length 8192. Only two 
roots of the Y transform of the Gaussian lie inside 
this radius and the value of the Y transform of the 
noise at these roots is quite small so that in the case 
of 24 bit accuracy even radiused FFT inversion of- 
fers a good estimate of the answer if adequate zero 
padding is provided. 



Figure 1 shows the range noise spectrum versus 
the spectrum of the range signal-plus-noise after 
root projection on the circle of radius 0.997. 

Since the projected output lies in the range of the 
Gaussian under convolution, FFT division on a 
space of dimension 800 or greater is now accurate 
and each of the FFT components can be considered 
as lying in an orthogonal channel.^ 

Using the principle that throwing out a channel 
throws out both the signal and the noise in that 
channel, we have a simple criterion that we should 
reject all channels where the signal/noise ratio is 
less than one or, crudely, where the (signal-plus- 
noise)/noise ratio is less than V2. In this case that 
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Figure 1. 800 point FFT on the circle of radius 0.997 of the range signal-plus-noise and noise for 
the convolution of the time series in figures 2.4a[l] and 1[2] with 24 bit roundoff noise added 
(Example la). 



^ The lack of independence of the FFT components in the image 
space means that projecting out some components will produce 
inverses no longer zero beyond point #700, thus no longer satis- 
fying the multiplication of FFT's rule of the space of dimension 
M. Signal-to-noise analysis, then, is not rigorously correct in 
RPD on a channel by channel basis. To avoid this complication, 
filtering for RPD will be carried out in the domain space after 
the inversion. This change of dimension before filtering implies 
a priori information about the carry-over of "relevant" frequency 
features between FFT representations of different cardinality. 
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means essentially every channel is acceptable so 
that no filtering of the inverse should be necessary. 
Figure 2 shows a comparison of the errors of the 
estimated inverse versus correct answer for RPD, 
radiused FFT inversion, and SVD inversion using a 
subspace of dimension 100. As can be seen, the 
RPD error is the smallest, merging with the FFT 
results for higher frequencies. The spectrum of the 
RPD error is, of course, exactly the quotient of the 
spectra of the range noise with that of the Gaus- 
sian kernel (padded to 800 points which was the 
FFT dimension used for RPD). The significantly 
larger error associated with the best Lanczos/SVD 
estimate arises from the use of only a 100 dimen- 
sional Lanczos subspace. The algorithm chosen 
uses the signal-plus-noise range vector as one of 
elements of the Lanczos subspace, but i)The range 
basis obtained may not be adequate to expand the 
noise or ii) The range basis may not separate the 
signal and the noise well, iii) The domain basis may 
not be adequate to expand the answer accurately, 
and iv) The matrix used to represent the convolu- 



tion matrix has to be truncated at the last row per- 
turbing it from its correct form on the whole space. 

In this case the error appears to stem equally 
from the last two of these causes. The relative stan- 
dard deviation (r.s.d.) from the answer caused by 
inadequate domain basis is 1.6% while the r.s.d. for 
the best SVD estimate is 2,2%, Tne r.s.d. in ex- 
panding the noise is negligible. The two errors 
stand in the ratio about V2 which would be ex- 
pected from equal independent contributions. By 
comparison the r.s.d. of the error of the RPD esti- 
mate is 0.16% while that of the FFT estimate is 
0.25%, The relatively large error in the Lanczos/ 
SVD estimate for this example remains no matter 
how small the roundoff error; only when the inver- 
sion error becomes small relative to the induced 
domain noise does the Lanczos/SVD method be- 
come effective. The SVD inversion estimate is 
shown in figure 3; the FFT and RPD estimates are 
visually identical to the correct answer. 

The ill-conditioned nature of the deconvolution 
process becomes apparent when the output has 
more noise. A comparison of range signal-plus- 
noise to noise after RPD is given in figure 4 for the 
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Figure 2. Comparison of spectra of errors for 3 deconvolution methods for Example la. An 800 
point FFT was used on the circle of radius 0.997. 
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Figure 3. Example la: Overplot of the SVD inversion estimate and the correct answer. 



case where the output has been truncated to 8 bit 
accuracy. The same figure made before root pro- 
jection would appear visually to be substantially the 
same. The noise level is lowered about 0.6 dB by 
RPD. These spectra were also carried out on the 
circle of radius 0.997. Because the time curves were 
not renormalized after radiusing, the average noise 
level reads slightly below the -64 dB expected 
level for unradiused noise. As can be seen the (sig- 
naUplus-noise/noise) ratio is clearly greater tiian 
VZ until about 0.065 of the frequency band, drop- 
ping to less than one beyond 0.1. 

Figure 5 shows the comparison of the errors for 
the estimated inverses from radiused ITT inver- 
sion (again padded to 8192 points) and inversion 
following RPD. In this case radiused FFT inversion 
produces incorrect answers, not only in the esti- 
mate for the inverse, but in the range upon recon- 
volving the "inverse" with the kernel. The values of 



the Y transform of the noise at the two roots lying 
inside the radius 0.997 have now greatly perturbed 
the radiused FFT division result. 

The errors for the inverse estimated by RPD are 
determined by the ratio of the 800 point FFT of the 
root projected noise to the 800 point ITT of the 
zero-padded Gaussian. One sees that root projec- 
tion has preserved the low-frequency features of 
the signal which are contained in the channels of 
high signal/noise ratio. 

Even without knowing the true noise distribution 
we can estimate the meaningful channels and build 
a time domain approximation from them. Figure 6 
shows the complete spectrum of the root projected 
inverse. We can use obvious a priori information 
that all channels beyond 0.08 of the frequency 
range, where the steep rise begins, are dominated 
by noise. We use an 801 point optimal filter with a 
pass band of 0.0 to 0.065 and stopband of 0.08 to 
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Figure 4. 800 point FFT on the circle of radius 0.997 of the range signal-plus-noise and noise 
after RPD for the convolution of the time series in figures 2.4a [1] and 1 [2] with 8 bit roundoff 
noise added (Example lb). 



1.0 with attenuation of 90 dB, convolve this with 
the complete RPD inverse and "trim the ends" to 
select out the 700 points from point 401 to point 
1100. Figure 7 gives the graph of that filtered esti- 
mate from RPD. It has an overall r.s.d. of 10.5%, 
mostly due to errors near the front and back. 
Throughout much of the first 400 points, the r.s.d. 
is less than 3%. 

A signal-to-noise analysis similar to that for RPD 
can be carried out for SVD inversion. In this case 
the channels in which one carries out the analysis 
are no longer associated with the familiar fre- 
quency eigenvectors but are built by the SVD pro- 
cess and are peculiar to the kernel series. However, 
one has the same number of range dimensions as 
domain dimensions which makes analysis of the fil- 
tering process simpler in SVD than in RPD. 

Figure 8 presents the actual signal and noise 
spectra for the Lanczos/SVD representation in the 
case of 8 bit range data accuracy. The signal gener- 
ally falls with decreasing singular values while the 
noise remains relatively constant at the expected 
value of - 64 dB except for a rise for small singular 



values (large dimension number). The fraction of 
the noise which can be represented by only 100 of 
the 800 range eigenvectors is (799.995/800)2 of the 
total noise. This highlights one of the features asso- 
ciated with the Lanczos/SVD method, notably the 
linkage between the range data and the SVD basis 
whose first element is the signal-plus-noise data 
vector. Since the SVD basis vectors have a special 
form, however, it can be difficult to obtain informa- 
tion on the noise statistics for detailed analysis. 

The plots of signal-plus-noise and inverse SVD 
spectra are given in figure 9. As can be seen by 
comparing figures 8 and 9, it is harder to find a 
best guess from the available data shown in figure 9 
than it was for root projection. The technique of 
truncating above a singular value, which leads to a 
residual closest to the expected noise residual, is 
completely in error here. It gives rise to more than 
70 terms in the expansion and produces an esti- 
mate full of noise with about 1000% r.s.d. from the 
correct answer. However, based on an expected 
signal-plus-noise cutoff value of about - 61 dB, a 
reasonable guess might be to include all compo- 
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Figure 5. Comparison of errors for root projection and radiused FFT inversion for Example lb. 
An 800 point FFT was used on the circle of radius 0.997. 



nents from dimension 1 through 59.'' The use of 
this "reasonable" guess produces the estimate 
shown in figure 10. The result is remarkably similar 
to that obtained by filtering of the root projection 
inverse. It has an overall r.s.d. of 10.4%. 



3.2 Example II 

The second example concerns deconvolution in 
the presence of incomplete range information. The 
kernel for this example is the experimental curve 
shown in figure 2.4a of reference [1] which was the 
answer in the previous example [1]. This time only 
400 output data points are given as shown in figure 
11. The error added was -40 dB random noise, 



'' One could also use a S/N cutoff at —61 dB; the resultant esti- 
mate in this case is visually identical to that chosen. However, 
we have found by experience that something like a Gibbs phe- 
nomena also occurs in SVD expansions so that including all 
components up to a given dimension number tends to give 
slightly less oscillations in the estimates than using a S/N cutoff 
A regularization type approach to this problem would be to 
use all channels for the inversion, but to replace each singular 
value, A, by A +€ where e is some appropriately chosen small 
constant. 



which produces an r.s.d. of about -41 dB. The 
noise is clearly visible in the figure. The answer, 
shown in figure 12, represents a possible trans- 
ducer response function. The output data then rep- 
resents the response of the transducer to the same 
stimulus that produced the kernel response of fig- 
ure 2.4a in reference [1] for the calibrating capac- 
itive transducer. 

As would be expected in a problem with a rising 
transient front end, the forward substitution ap- 
proach is inapplicable, producing relative errors of 
5000 on the third point and overflowing soon there- 
after. The result of applying FFT division using an 
800 point FFT is shown in figure 13. The range 
r.s.d. for this curve is 158% indicating a completely 
erroneous result even though the correct answer 
has less than 200 non-zero points and 400 points 
were given. If zero padding is added to the FFT 
(8192 point FFT) the range r.s.d. drops to 44% 
rather than the 0.8% needed to be within the ex- 
pected range noise error. The domain r.s.d. against 
which the SVD and RPD estimates will be com- 
pared below is 1976%. Noting that only two roots of 
the kernel lie below 0.979, one can apply radiused 
FFT division with a radius of 0.975 to produce a 
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Figure 6. Spectrum from 800 pt FFT (with radius = 1.0) of estimated inverse for projection show- 
ing evident break between signal and noise dominated regions. 



range r.s.d. of 1.1%. The range r.s.d. is almost 
within acceptable limits, but the domain errors are 
substantially greater than those produced by either 
Lanczos/SVD or RPD. 

To deal with this model problem for which one 
has only limited a prion information (mostly that 
the answer should not be filled with high frequency 
noise), we introduce the cross-cutting deconvolu- 
tion algorithm (CCD) shown schematically in fig- 
ure 14. This algorithm uses the Lanczos/SVD and 
RPD methods in tandem to produce a more robust 
estimate which is at least as good as the particular 
choice of either estimate and often better than 
both of these methods. Lanczos/SVD and RPD 
with optimal filtering are first applied indepen- 
dently to produce first estimates to the decon- 
volved inverse. Typically, with a kernel whose 
spectrum is dominated by low frequency compo- 
nents, the SVD estimate will show some of the 
most prominent high frequency features, but will 
have reduced low frequency fidelity. 

The RPD estimate, on the other hand, will tend 
to have good low frequency fidelity but reduced 
resolution of some of the high frequency features 
and possibly greater end noise. The outputs from 



each of these algorithms, conservatively filtered to 
avoid extraneous features, is then fed into the 
other algorithm as an a priori estimate. Each of 
these two methods uses a different orthogonal de- 
composition to separate signal from noise and the 
reuse of the output from one of these algorithms 
provides essentially no improvement in the estimate 
(although RPD may be rerun to reduce end noise). 
However, what one algorithm may discard as noise 
can contain useful signal when decomposed using 
the other algorithm. The second, cross-cut esti- 
mates are then frequently significant improvements 
upon the initial estimates. They each, of course, 
contain induced domain noise which may be par- 
tially independent so that the average of the esti- 
mates can be expected to provide some overall 
improvement over the second approximations. 

Figure 15a shows the signal-plus-noise and do- 
main spectra and figure 15b the actual signal and 
noise spectra upon applying Lanczos/SVD inver- 
sion to the model problem data. Based on figure 
15a with an approximate signal-plus-noise value of 
-37 dB, reasonable choices might lie between se- 
lecting all components from 49 up to 80. As one 
increases the dimension number past 49, the actual 
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Figure 7. Filtered root projection estimate for Example lb. 



r.s.d. decreases more or less steadily (cf. fig. 15b) 
from 103% to 65% at dimension number 77. How- 
ever, as the domain noise increases (and due also 
to increasing "high frequency" characteristics of 
the SVD base vectors) more extraneous features 
become added to the signal. In this case because 
the eigenvalues are not too small in the 49 to 80 
range vis-a-vis those below 49, the "risk" is not too 
great of introducing a great deal of "high fre- 
quency" domain noise by inadvertently adding a 
channel with excess noise. The only guideline then 
for choosing among these candidate estimates is a 
priori information about the signal.^ We choose the 



^The Lanczos dimension used throughout this paper is 100. 
However, a special run was made for this example using Lanc- 
zos dimension 200. Although the error in representation of the 
answer was much better in the 200 case, the range noise was 
sufficiently large that only small improvement occurred in the 
r.s.d. of the estimates. The main improvement was in the signal- 
plus-noise diagram which showed a distinct plateau now making 
the choice of dimension number 77 as a cutoff almost apparent. 



estimate made up from all dimension numbers 
through 62. This estimate has an overall r.s.d. of 
75%. No matter which of these estimates is used as 
a starting point, the final output of the cross-cut- 
ting algorithm varies only between 52% and 54% 
r.s.d. 

The RPD techniques for finding the initial ap- 
proximation from incomplete range data differ 
somewhat from those of the SVD method. The po- 
tentially large amount of noise that can be intro- 
duced by the a priori part of the range estimate 
required by the RPD method can preclude a sig- 
nal-to-noise analysis similar to that used for SVD. 
Instead we examine the a posteriori data separately 
to estimate bounds on the useful part of the fre- 
quency range and re-run RPD on windowed raw 
estimates to remove excessive end noise. 

Figure 16a shows the FFT of the 400 points of a 
posteriori data after they have been windowed using 
a 400 point window built from a maximal ripple 
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Figure 8. The true signal and noise spectra for Example lb using the SVD basis (not the Fourier 
basis) on a Lanczos subspace of dimension 100. 



filter to smooth the frequency distribution. The fig- 
ure shows essentially no information above the 
- 40 dB noise level beyond 30% of the Nyquist fre- 
quency. Any information about the inverse above 
the relative frequency of 0.3, then, must either be a 
priori information or information obtained from the 
SVD estimate. In other cases there will be informa- 
tion over the noise level through high frequency or 
other "difficult" parts of the frequency range 
where the kernel has a small spectral amplitude. 

With no a priori information available we employ 
a simple extrapolation procedure to create the in- 
put range signal for the first RPD approximation. 
Figure 16b shows the normalized signal-plus-noise 
and domain spectra after RPD inversion of the 
given data with the endpoint values continued out 
to the end of the range data and then windowed. A 
time weighting of 1:16 (one in the given data and 
16 in the extrapolated data) was used to bring the 
range residuals to well within the expected -41 dB 
r.s.d. over the initial data. When the weight factor 
is increased much beyond 20 in this example, expo- 
nentially increasing noise terms begin to show up in 
the inverse and negligible further information 
seems to be obtainable. In other less noisy or ill- 



conditioned cases the weight factor can be in- 
creased to several thousand. 

Figure 16c shows the spectra of the true signal 
and range noise for the first RPD approximation. 
The range signal is merely the product of the 800 
point FFT's of the kernel and the answer. The 
range noise, on the other hand, is not simply the 
- 40 dB fl posteriori noise, but is here completely 
dominated by the a priori noise which, because of 
the smooth extrapolation used, has very large low 
frequency components. Without a priori informa- 
tion, there seems to be no simple way of estimating 
the actual range noise which can be dominated by 
such errors. The initial procedure for RPD in this 
case, then, is to choose as broad-band a filter as 
possible producing an estimate which doesn't show 
excessive high frequency noise and to use the gen- 
eral a priori information that the inverse should not 
have a rising high frequency spectrum. No matter 
which filter is chosen, however, one can expect ex- 
tensive low frequency errors. How they manifest 
themselves is shown in figure 17. 

Figure 17 shows the SVD and RPD first esti- 
mates. As described above the overall r.s.d. for the 
SVD estimate is 75%. The overall r.s.d. for the 



358 



Volume 96, Number 3, May-June 1991 

•Journal of Research of the National Institute of Standards and Technology 



UO.OO' 



0.00 



tn 



CD 



o -40. 00' 



o 



a. 

CE 



-80. 00' 



120.00 



160.00 



J 



^ 




A / ■. !■ ■ 







LEGEND 
RfiNGE SlGNfiL + NOISE (SVG) 
DOMRIN INVERSE (SVD) 



0.00 



20.00 40.00 60.00 80.00 
LfiNCZOS DIMENSION NUMBER 



100.00 



Figure 9. Range signal-plus-noise and estimated inverse spectra for Example lb using the appro- 
priate SVD bases built in the Lanczos subspace of dimension 100. 



RPD estimate is 471%. Examination of the two 
curves shows this difference to be due to the much 
greater end error in the RPD estimate. Both of 
these methods tend to produce larger errors at the 
end of the time domain since the end values of the 
answer are involved in fewer terms of the a posteri- 
ori data. RPD, particularly without a reasonable a 
priori range guess, produces such errors due to 
leakage of the large a priori errors into the a poste- 
riori part of the time domain.* It's easy to confirm 
in this case that there is almost no valid data be- 
yond point 140 in either estimate. This is done for 
RPD by using the RPD estimate of figure 17 win- 
dowed to 140 points as an a priori estimate and 
reinverting. A lower time weighting ratio of 1:4 is 
all that's needed in this case to produce good resid- 
uals since the a priori range data are more accu- 
rate. After filtering the resultant estimate is almost 



" Errors also tend to occur very near the causal time boundary in 
both methods when extensive projection is employed. One can 
give a heuristic interpretation of this front end error for RPD 
based on the fact that the RPD estimate is low-pass filtered. By 
reducing the frequency range by a factor of Af, say, lower reso- 
lution occurs in the time domain which shows up most clearly in 
the first M points of the time domain. 



indistinguishable from the original for the first 140 
points and has only small oscillations past point 
140. This indicates that the large end oscillations in 
the initial estimate were forced by a priori range 
errors. The overall r.s.d. of the error of this RPD 
estimate is 60%. 

A similar process to show that there is almost no 
valid data beyond point 140 of the SVD estimate 
can also be carried out. However, there is nothing 
intrinsically "wrong" with the estimate produced 
by the RPD process which returns range residuals 
almost within the expected noise limits (some in- 
crease in residuals is due to filtering). It is only that 
the a priori range data chosen forced features in 
the answer which were irrelevant to the essential 
features required by the a posteriori data. In fact 
one can establish a different "model" problem in 
which the correct answer is the curve given as the 
RPD estimate in figure 17. Application of SVD to 
400 points of range data from this second "model" 
problem with similar noise added produces esti- 
mates similar to those of the model problem as we 
have presented it with only slightly increased end 
oscillations. The overall r.s.d. of the error of the 
similar 62 component estimate from the answer to 
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Figure 10. Best Lanczos/SVD estimate for Example lb using S/N projection. 



the alternate problem, for instance, is 366%. The 
extra data needed to identify the end features of 
the alternate answer are contained in the last 399 
points of missing data. 

In most problems with 50% incomplete range 
data one can still obtain accurate information 
about the answer through most of the time domain 
involved. Typical examples we have checked show 
good accuracy for 60% to 95% of the time range. 
An assessment of the importance of missing end 
data may be possible in RPD using knowledge of 
the positions of the roots in the root transform of 
the answer vis-a-vis those of the transform of the 
kernel. 

The second stage of the cross-cut algorithm uses 
each of the two first stage outputs as a priori inputs 
to the alternate algorithm. For this example, since 
we have established that little meaningful data lie 
beyond point 140, we will use the windowed first 
estimates rather than the full estimates as a priori 



inputs. For the SVD process, then, the windowed 
RPD first estimate is convolved with the kernel 
and the first 400 points of the resultant time series 
are subtracted from the a posteriori data. Lanczos/ 
SVD deconvolution is then performed on this re- 
duced problem and the estimated inverse added 
back onto the a priori estimate to produce the sec- 
ond estimate. 

Figure 18a shows the signal-plus-noise and in- 
verse SVD spectra and figure 18b shows the true 
signal and noise spectra for this reduced SVD 
problem. With a - 37 dB signal-plus-noise approxi- 
mate cutoff, a reasonable choice from figure 18a 
would be a cutoff at 12 vectors, which in this case is 
confirmed by the curves in figure 18b. The r.s.d. for 
the second SVD estimate is then 56%. 

Figure 19a shows the signal-plus-noise and in- 
verse RPD spectra and figure 19b shows the actual 
range signal and noise spectra for the reduced 
RPD problem which has as its a priori input the 
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Figure 11. Synthetic output data with noise for model trans- 
ducer calibration problem (Example II). 

windowed first estimate from the Lanczos/SVD al- 
gorithm. The extensive high frequency a priori 
noise level is no longer present, and the apparent 
cutoff from figures 19a and 19b is between 0.16 and 
0.2 relative frequency. We choose the same filter 
(transition between 0.16 and 0.3 relative fre- 
quency) for this stage as in the first stage, since no 
high frequency noise is apparent. This second RPD 
estimate has an overall r.s.d. of 60%. 

The final step in the cross-cut algorithm is to av- 
erage the two estimates arising from the second ap- 
proximations. The resultant average is shown 
compared with the correct answer in figure 20. It 
has an overall r.s.d. of 52%. Since the geometric 
mean of the answer is 0.014, the actual domain 
noise is -42.5 dB or 2.5 dB less than the added 
range noise. Figure 21a gives a comparison of the 
spectrum of the cross-cut estimate, windowed to 
200 points to smooth the spectrum slightly, with 
that of the answer. Figure 21b gives a similar com- 
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Figure 12. Correct answer for Example II. 

parison of the phases. As can be seen from these 
figures, almost all the error in the estimate is at 
frequencies above 0.1. 



4. Summary 

A comparative development has been given of a 
number of related deconvolution methods includ- 
ing a new algorithm based upon root projection. 
Radiused FFT division, singular valued decomposi- 
tion (SVD) and root projection deconvolution 
(RPD) were compared in detail in two examples. 
The first example involved deconvolution of a 101 
point Gaussian filter from a relatively complex 800 
point waveform in the presence of noise. Two lev- 
els of noise were used. The concept of signal-to- 
noise filtering was introduced and applied to both 
SVD and RPD. It was shown to be particularly 
powerful when used in conjunction with the SVD 
method. In both cases RPD provided the most ac- 
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Figure 13. Result of using 800 point FFT division to estimate 
answer for Example II. 

curate inversion estimate although this accuracy 
was matched by the SVD method for the higher 
noise level. The second example was a model trans- 
ducer calibration problem with typically incomplete 
range data. A robust cross-cut deconvolution 
(CCD) algorithm utilizing the orthogonal decompo- 
sitions of both the RPD and SVD methods was used 
in this example. The deconvolved estimate showed 
virtually no errors for frequencies below 10% of the 
Nyquist frequency. 

Root transforms were applied in this work only 
for one dimensional time series. The study of the 
comparable Riemann surfaces for higher dimen- 
sional series or the use of boundary roots in finding 
signatures for the Radon transforms of such series 
have yet to be studied. 
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Figure ISa, Range signal-plus-noise and estimated inverse spectrum for Example 
II using the appropriate SVD bases (not the Fourier bases) built in the Lanczos 
subspace of dimension 100. 
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Figure ISb. The true signal and noise spectra for Example II using the SVD range 
basis on the Lanczos subspace of dimension 100. 
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Figure 16a. Normalized 400 point FFT of the windowed output data of figure 11. 
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Figure 16b. Signal-plus-noise and estimated inverse spectra for example II using an 
800 point FFT on the result of root projection where the data of figure 11 were 
simply end extrapolated out to 800 points. 
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Figure 16c. Spectra of 800 point FFT of true range signal for Example II and noise 
from initial RPD inversion. 
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Figure 17. Initial SVD and RPD estimates for Example II. 
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Figure 18a. Range signal-plus-noise and estimated inverse spectra for the second 
stage of the cross-cut algorithm in Example II. The appropriate SVD bases built in 
the Lanczos subspace of dimension 100 are used. The input is the first stage RPD 
estimate windowed to 140 points. 
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Figure 18b. The true signal and noise spectra for the second stage of the cross-cut 
algorithm in Example II using the SVD range basis on the Lanczos subspace of 
dimension 100. 
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Figure 19a. Signal-plus-noise and estimated inverse spectra for the second stage of 
the cross-cut algorithm in Example II. An 800 point FFT is used on the result of 
root projection from the windowed first stage SVD input. 
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Figure 19b. The true range signal and noise spectra for the second stage of the 
cross-cut algorithm in Example II. An 800 point FFT is used on the result of root 
projection from the first stage SVD input. 
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Figure 20. The CCD estimate for Example II. 
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Figure 21a. Spectrum of the CCD estimate for Example II. 
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Figure 21b. Phase of the CCD estimate for Example II. 
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